 #This script performs OLS and quantile (linear) regressions for 4 species of skinks, and prepares a datafile to compare residual abundance to microclimate

#Clear workspace

rm(list=ls())

in.dir = "C:/R/In/"
out.dir = "C:/R/Out/Georef Scatterplots/avgabund/"

#Set working directory
setwd(in.dir)

#load necessary libraries

library("SDMTools")
library("quantreg")
library("adehabitat")

#Import the dataset for analysis

samples= read.csv("samplesmaster.csv",header=T)
georef = read.csv("georefmaster.csv", header=T)
s.fields = names(samples)
g.fields = names(georef)

#First begin by working with the samples data, plotting relationships between all variables and sample abundance
#Use a for loop to perform both quantile and OLS regression between all variables and samples abundance

num=1

for(fields in g.fields)

  {
  
  t.data=data.frame(georef$cr_avgabund, georef[,num])

  png(paste(out.dir,"CARRUBR","/",fields,".png", sep=""), height=800, width=800)

  plot(t.data[,2], t.data[,1], xlab=fields, ylab="avgabund", main="CARRUBR")
  
  t.rq1 = rq(t.data[,1]~t.data[,2], tau=.90, data=t.data, method="br", model=T)
  t.rq2 = rq(t.data[,1]~t.data[,2], tau=.95, data=t.data, method="br", model=T)
  t.rq3 = rq(t.data[,1]~t.data[,2], tau=.975, data=t.data, method="br", model=T)
  t.rq4 = rq(t.data[,1]~t.data[,2], tau=.99, data=t.data, method="br", model=T)
  t.ols = lm(t.data[,1]~t.data[,2], data=t.data)

  abline(t.rq1, col="red")
  abline(t.rq2, col="red")
  abline(t.rq3, col="red")
  abline(t.rq4, col="red")
  abline(t.ols, col="blue")

  
  dev.off()
  
  num=num+1

  }
  
num=1

for(fields in g.fields)

  {

  t.data=data.frame(georef$gq_avgabund, georef[,num])

  png(paste(out.dir,"GNYQUEE","/",fields,".png", sep=""), height=800, width=800)

  plot(t.data[,2], t.data[,1], xlab=fields, ylab="avgabund", main="GNYQUEE")

  t.rq1 = rq(t.data[,1]~t.data[,2], tau=.90, data=t.data, method="br", model=T)
  t.rq2 = rq(t.data[,1]~t.data[,2], tau=.95, data=t.data, method="br", model=T)
  t.rq3 = rq(t.data[,1]~t.data[,2], tau=.975, data=t.data, method="br", model=T)
  t.rq4 = rq(t.data[,1]~t.data[,2], tau=.99, data=t.data, method="br", model=T)
  t.ols = lm(t.data[,1]~t.data[,2], data=t.data)

  abline(t.rq1, col="red")
  abline(t.rq2, col="red")
  abline(t.rq3, col="red")
  abline(t.rq4, col="red")
  abline(t.ols, col="blue")


  dev.off()

  num=num+1

  }

num=1 
for(fields in g.fields)

  {

  t.data=data.frame(georef$sb_avgabund, georef[,num])

  png(paste(out.dir,"SAPBASI","/",fields,".png", sep=""), height=800, width=800)

  plot(t.data[,2], t.data[,1], xlab=fields, ylab="avgabund", main="SAPBASI")

  t.rq1 = rq(t.data[,1]~t.data[,2], tau=.90, data=t.data, method="br", model=T)
  t.rq2 = rq(t.data[,1]~t.data[,2], tau=.95, data=t.data, method="br", model=T)
  t.rq3 = rq(t.data[,1]~t.data[,2], tau=.975, data=t.data, method="br", model=T)
  t.rq4 = rq(t.data[,1]~t.data[,2], tau=.99, data=t.data, method="br", model=T)
  t.ols = lm(t.data[,1]~t.data[,2], data=t.data)

  abline(t.rq1, col="red")
  abline(t.rq2, col="red")
  abline(t.rq3, col="red")
  abline(t.rq4, col="red")
  abline(t.ols, col="blue")


  dev.off()

  num=num+1

  }

num=1  
for(fields in g.fields)

  {

  t.data=data.frame(georef$lc_avgabund, georef[,num])

  png(paste(out.dir,"LAMCOGG","/",fields,".png", sep=""), height=800, width=800)

  plot(t.data[,2], t.data[,1], xlab=fields, ylab="avgabund", main="LAMCOGG")

  t.rq1 = rq(t.data[,1]~t.data[,2], tau=.90, data=t.data, method="br", model=T)
  t.rq2 = rq(t.data[,1]~t.data[,2], tau=.95, data=t.data, method="br", model=T)
  t.rq3 = rq(t.data[,1]~t.data[,2], tau=.975, data=t.data, method="br", model=T)
  t.rq4 = rq(t.data[,1]~t.data[,2], tau=.99, data=t.data, method="br", model=T)
  t.ols = lm(t.data[,1]~t.data[,2], data=t.data)

  abline(t.rq1, col="red")
  abline(t.rq2, col="red")
  abline(t.rq3, col="red")
  abline(t.rq4, col="red")
  abline(t.ols, col="blue")


  dev.off()

  num=num+1

  }